subroutine PatchTestQ8
    use Mesh2DQ8_m
    use Mesh2DQ4_m
    use FES2DQ82_m
    use SMatrix_m
    use IterSSolver_m
    implicit none
    type(mesh2DQ4) :: mesh2D0
    type(mesh2DQ8) :: mesh2D
    type(FES2DQ82) :: fs
    type(SMatrix) :: K
    type(IterSSolver) :: solver
    real*8, allocatable :: x(:), RHS(:), u(:), v(:)
    real*8, allocatable :: ex(:), ey(:), exy(:), tx(:), ty(:), txy(:), fx(:), fy(:)
    
    real*8, parameter :: E = 1.0D6, Nu = 0.25D0
    integer :: i
    integer, parameter :: BDLabel(1) = (/4/)
    character*80, parameter :: dir = "PatchTest\"
    character*80, parameter :: meshFile = trim(dir)//"PatchTestMeshQ8.mesh"
    character*80, parameter :: outputFile = trim(dir)//"solutionQ8.dat"
    
    !call mesh2D.readMesh_msh(meshFile)
    !call printNodesQ8(mesh2D)
    !call squareMeshQ8(mesh2D, 5, 5)
    !call exportMesh(outputFile, mesh2D)
    !call mesh2D0.square(5, 5)
    !call mesh2D.initMesh(mesh2D0)
    call mesh2D.square(5, 5)
    !initMesh(mesh2D0)
    call fs.initFES(mesh2D)
    call K.initSMatrix(fs%nDof)
    allocate(x(fs%nDof), RHS(fs%nDof))
    x = 0.0D0; RHS = 0.0D0
    allocate(u(fs%nDof / 2), v(fs%nDof / 2))
    allocate(ex(fs%nDof / 2), ey(fs%nDof / 2), exy(fs%nDof / 2), tx(fs%nDof / 2), ty(fs%nDof / 2), txy(fs%nDof / 2), fx(fs%nDof / 2), fy(fs%nDof / 2))
    
    call fs.assembleElasticityStiff(E, Nu, 0, K)
    !call K.sort()
    !call K.display()
    call fs.assembleSurfaceForce(100.0D0, 0.0D0, 2, 0, RHS)

    call fs.applyBC(K, RHS, BDLabel, 0.0D0, 0.0D0, 1)
    call solver.initIterSSolver()
    call solver.PCGSSOR(K, RHS, x)
    call fs.extractDisplacement(x, u, v)
    call fs.computeStrainStressInternalForce(u, v, E, Nu, ex, ey, exy, tx, ty, txy, fx, fy, 0)
    
    !call exportMeshSolution(outputFile, mesh2D, u, v)
!    write(*, "(F10.6)") x
!    deallocate(x, RHS)
    
    write(*, *) "Displacement: "
    write(*, *) "u min:", minval(u), "max:", maxval(u)
    write(*, *) "v min:", minval(v), "max:", maxval(v)

    write(*, *) "Strain: "
    write(*, *) "ex min:", minval(ex), "max:", maxval(ex)
    write(*, *) "ey min:", minval(ey), "max:", maxval(ey)
    write(*, *) "exy min:", minval(exy), "max:", maxval(exy)

    write(*, *) "Stress: "
    write(*, *) "tx min:", minval(tx), "max:", maxval(tx)
    write(*, *) "ty min:", minval(ty), "max:", maxval(ty)
    write(*, *) "txy min:", minval(txy), "max:", maxval(txy)

!    write(*, *) "d[x]       d[y]       |  fint[x]    fint[y]    |  sigma[xx]   sigma[yy]   tau[xy]"
!    do i = 1, fs%nDof / 2
!        write(*, "(7E13.3)") u(i), v(i), fx(i), fy(i), tx(i), ty(i), txy(i)
!    end do
    
    call K.deallocateSMatrix()
    call mesh2D.deallocateMesh()
contains
real*8 function uFunc(x, y, time)
    implicit none
    real*8, intent(in) :: x, y, time
    uFunc = (x + 0.5 * y) *0.001D0
end function

real*8 function vFunc(x, y, time)
    implicit none
    real*8, intent(in) :: x, y, time
    vFunc = (y + 0.5 * x) *0.001D0
end function

end subroutine